Protocol to get the hydraulic conductivities along the root axes of maize root.

The results of coupling GRANAR with MECHA is the estimation of the hydraulic conductivities of a root. GRANAR will generate root cross section anatomy. MECHA will use the synthetic root cross section to calculate the radial conductivities (Kr) and the axial conductance (Kx).

To access the information needed to run GRANAR, it requires free-hand section or permanent root cross section. Staining the cross section will allow to spot the different apoplastic barriers.

In this example, free-hand section were stain with aqueouse solution of berberin (1h) and post-stain with analine blue (30min).

Basal root at 15cm from the tip stain with berberin and post-stain with analine blue

Basal root at 15cm from the tip stain with berberin and post-stain with analine blue

# Dependencies

library(tidyverse)
library(cowplot)
library(readxl)
library(plyr)
library(deldir)
library(alphahull)
library(sp)
library(xml2)
library(GGally)
# List of function from the GRANAR packages that have been slighly modified and functions for this script 
source("R/granar.R")
`%!in%` <- compose(`!`, `%in%`)

Loading experimental data

root_name <- c("SBR", "Tap", "Basal", "LAT", "LAT_A")

Data_params <- read_excel("Data_params.xlsx")%>%
  as.tibble()
Data_params$cell[Data_params$Root_id == "05_01" & Data_params$tissue == "xylem" & !is.na(Data_params$cell)] <- Data_params$cell[Data_params$Root_id == "05_01" & Data_params$tissue == "xylem" & !is.na(Data_params$cell)]*2
params <- read_param_xml("www/Zea_mays_B73_S_01.xml")

# number of time the generation of anatomies will be replicated.
n_rep <- 2 # always >1

The anatomical features ploted against the distance to the tip

We observed that many of them follow a linear regression.

Data_params%>%
  filter(tissue %!in% c("proto", "xylem"),
         !is.na(tissue),
         tissue != "aer")%>%
  ggplot(aes(dist, cell))+
  geom_point(alpha = 0.2, aes(colour = Root))+
  geom_smooth(method = "lm", aes(colour = Root))+
  facet_wrap(~tissue)+
  ylab("Cell size [um]")+
  xlab("distance from the tip [cm]")+
  viridis::scale_colour_viridis(discrete = T)

Data_params%>%
  filter(tissue %in% c("cortex"))%>%
  ggplot(aes(Root, layer))+
  geom_boxplot(aes(colour = Root))+
  facet_grid(~tissue)+
  ylab("layer width [um]")+
  xlab("distance from the tip [cm]")+
  viridis::scale_colour_viridis(discrete = T)

Data_params%>%
  filter(tissue == "xylem",
         !is.na(tissue))%>%
  ggplot(aes(dist, cell))+
  geom_point(alpha = 0.2, aes(colour = Root))+
  geom_smooth(method = "lm", aes(colour = Root))+
  facet_wrap(~tissue)+
  ylab("Cell size [um]")+
  xlab("distance from the tip [cm]")+
  viridis::scale_colour_viridis(discrete = T)

stelar <- Data_params%>%
  filter(tissue %in% c("stele"))

stelar%>%
  ggplot(aes(dist, layer))+
  geom_point(aes(colour = Root))+
  geom_smooth(aes(colour = Root), method = "lm")+
  facet_grid(~tissue)+
  viridis::scale_colour_viridis(discrete= T)

xylem <- Data_params%>%
  filter(tissue %in% c("xylem"))

xylem %>%
  ggplot(aes(dist, cell, colour = Root))+
  geom_point()+
  geom_smooth(method = "lm")+
  viridis::scale_colour_viridis(discrete= T)

xylem %>%
  ggplot(aes(dist, layer, colour = Root))+
  geom_point()+
  geom_smooth(method = "lm")+
  viridis::scale_colour_viridis(discrete= T)

cortex <- Data_params%>%
  filter(tissue %in% c("cortex"))

cortex %>%
  ggplot(aes(dist, layer, colour = Root))+
  geom_point()+
  geom_smooth(method = "lm")+
  ylab("cortex width")+
  viridis::scale_colour_viridis(discrete= T)

The correlation between xylem and stele

cor_ma_tmp <- Data_params%>%
  filter(tissue %in% c("xylem"))%>%
  dplyr::group_by(Image)%>%
  dplyr::summarise(xylem_size = mean(cell))%>%
  ungroup()
cor_ms_tmp <- Data_params%>%
  filter(tissue %in% c("stele"))%>%
  dplyr::group_by(Image)%>%
  dplyr::summarise(stele_size = mean(layer))%>%
  ungroup()

cor <- merge(cor_ma_tmp, cor_ms_tmp, by = "Image")

ggpairs(cor%>%select(-Image))

The location of the apoplastic barriers formation and the maturation of the metaxylem.

rom <- read_excel("Root_maturation.xlsx")%>%
  as.tibble()

Data_params%>%
  filter(!is.na(Xyl_m) )%>%
  ggplot()+
  geom_point(aes(dist, Xyl_m, colour = Root))+
  geom_line(aes(dist, Xyl, colour = Root), size = 1, alpha = 0.8, data = rom)+
  facet_wrap(~Root, scale = "free")+
  viridis::scale_colour_viridis(discrete= T)

Data_params%>%
  filter(!is.na(Xyl_m), Root_id != "14_01" )%>% # outliers
  ggplot()+
  geom_point(aes(dist, -Apo, colour = Root))+
  geom_line(aes(dist, -Apo, colour = Root), size = 1, alpha = 0.5, data = rom%>%filter(!is.na(Apo)))+
  facet_wrap(~Root, scale = "free")+
  viridis::scale_colour_viridis(discrete= T)

Gathering data for an overall view of the root cross sections.

# Prepare data
Data <- Data_params%>%
  dplyr::group_by(Image, tissue, dist, Root, Root_id)%>%
  dplyr::summarise(cell = mean(cell, na.rm = T),
                   layer = mean(layer, na.rm = T))%>%
  ungroup()%>%
  filter(!is.na(tissue))  

stele = CT = nX = xylem = CF <- NULL

for (im in unique(Data$Image)) {
  # Get the total diameter of the root
  
  tmp <- Data%>%
    filter(Image == im)
  if(nrow(tmp) == 8){
    
    tot <- tmp$cell[tmp$tissue == "endo"]+
      tmp$cell[tmp$tissue == "pericycle"]+
      tmp$cell[tmp$tissue == "exo"]+
      tmp$cell[tmp$tissue == "epi"]+
      tmp$layer[tmp$tissue == "stele"]/2+
      tmp$layer[tmp$tissue == "cortex"]
    
    total <- tibble(Image = im, 
                    tissue = c("total"),
                    cell = NA, 
                    layer = tot, 
                    dist = unique(tmp$dist), 
                    Root = unique(tmp$Root),
                    Root_id = unique(tmp$Root_id))
    stele <- cbind(stele, tmp$layer[tmp$tissue == "stele"]/2+tmp$cell[tmp$tissue == "endo"]+
      tmp$cell[tmp$tissue == "pericycle"])
    CF <- cbind(CF, tmp$layer[tmp$tissue == "cortex"]/tmp$cell[tmp$tissue == "cortex"])
    CT <- cbind(CT, tot)
    nX <- cbind(nX, tmp$layer[tmp$tissue == "xylem"] )
    xylem <- cbind(xylem, tmp$cell[tmp$tissue == "xylem"])
    
    Data <- rbind(Data, total)
  }
}
sct <- tibble(stele= c(stele), CT= c(CT), nX = c(nX), xylem = c(xylem), CF = c(CF))

correlation between xylem area and stele area under log.

## 
## Call:
## lm(formula = ln_XSA ~ ln_SXA, data = sct)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.62525 -0.09282  0.00842  0.09605  0.43372 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.34534    0.12325  -19.03   <2e-16 ***
## ln_SXA       1.22717    0.01344   91.29   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1649 on 95 degrees of freedom
## Multiple R-squared:  0.9887, Adjusted R-squared:  0.9886 
## F-statistic:  8334 on 1 and 95 DF,  p-value: < 2.2e-16

GRANAR

The anatomical data can be modelled to set the GRANAR parameter. The function Get_granarparam take the anatomical data for each root type individualy and look for a relation between the feature and the distance from the root apex. if p-value > 0.05, then the average value for the feature is took.

simple <- NULL
mapp <- c(seq(0.5, 2, by = 0.25),3:15,20,30,40)
# root_name <- root_name[-1]
for (root_type in root_name){
  mapp <- c(seq(0.5, 2, by = 0.25),3:15,20,30,40)
  # Shoot Born root parameter
  param <- Get_granarparam(Data_params = Data_params%>%
                             filter(Root == root_type,
                                    tissue != "aer"))
  tiss_param <- param[[1]]
  layer_param <- param[[2]]

  # Run GRANAR on 30 cm and make 3 repetitions
  all_roots = nodes <- NULL
  for (d in mapp){
    message(paste0(">>>>> Generation of ",root_type," cross section at ",d,"cm from tip. <<<<<"))

    for(rep in c(1:n_rep)){
      print(rep)

      params <- set_params(tiss_param, layer_param, params)

      # n metaxylem vessels
      nX <- layer_param$one[layer_param$lay =="xylem"]+d*layer_param$slope[layer_param$lay =="xylem"]
      SXA <- pi*((layer_param$one[layer_param$lay =="stele"]+d*layer_param$slope[layer_param$lay =="stele"])/2)^2
      XSA <- exp(xyl_int+xyl_m*log(SXA,base = exp(1)))
      xyl_area <- XSA/nX
      # xylem size
      c_xylem <- sqrt(xyl_area/pi)
      c_stele <- tiss_param$one[tiss_param$tiss=="stele"]+d*tiss_param$slope[tiss_param$tiss =="stele"]
      params$value[params$name == "xylem" & params$type == "max_size"] <- c_xylem/1000-c_stele/1000

      sim <- create_anatomy(parameters = params)
      if (d %in% c(2, 5, 10 )){print(plot_anatomy(sim))}

      if(file.exists("MECHA_GRANAR/cellsetdata/B73/current_root.xml")){
        file.remove("MECHA_GRANAR/cellsetdata/B73/current_root.xml")
      }
      write_anatomy_xml(sim = sim, path = "MECHA_GRANAR/cellsetdata/B73/current_root.xml")

      fc <- file.copy(from = "MECHA_GRANAR/cellsetdata/B73/current_root.xml",
                        to = paste0("MECHA_GRANAR/cellsetdata/B73/",root_type,"_", d,"_",rep,".xml"),
                        overwrite = T)
      all_roots <- rbind(all_roots, sim$output%>%
                             mutate(dist = d,
                                    repet = rep))
      nodes <- rbind(nodes, sim$nodes%>%
                         mutate(dist = d,
                                repet = rep))
      }
  }
    #Save data
  write.csv(all_roots, paste0("all_roots_",root_type,".csv"))
  write.csv(nodes, paste0("nodes_",root_type,".csv"))
  
  simple <- rbind(simple, 
                  tiss_param%>%transmute (intercept = one,
                                          slope = slope,
                                          tissue = tiss,
                                          type = "cell diameter",
                                          root = root_type),
                  layer_param%>%transmute (intercept = one,
                                          slope = slope,
                                          tissue = lay,
                                          type = "layer feature",
                                          root = root_type))
  
}

write.csv(simple, "linear_params.csv")

Summary of the GRANAR runs over the moddeled root anatomical features.

nodes%>%
  filter(repet == 1,
         dist %in% c(5,10,15,20),
         name %in% c("Tap", "SBR", "Basal"))%>%
  ggplot()+
  geom_segment(aes(x = x1, xend = x2, y = y1, yend = y2))+
  theme_classic()+
  coord_fixed()+
  facet_grid(name~dist)

nodes%>%
  filter(repet == 1,
         dist %in% c(1,2,5,7,10),
         name %in% c("LAT", "LAT_A"))%>%
  ggplot()+
  geom_segment(aes(x = x1, xend = x2, y = y1, yend = y2))+
  theme_classic()+
  coord_fixed()+
  facet_grid(name~dist)+
  ylim(0,0.2)

Benchmark

Does the simulated cross section have similar anatomical feature as the one that were used to set the parameters.

diameter <- nodes%>%
  dplyr::group_by(name, dist, repet, type)%>%
  dplyr::summarise(min_x = min(x1),
                   min_y = min(y1),
                   max_x = max(x1),
                   max_y = max(y1),
                   d = ((max_x-min_x)+(max_y-min_y))/2)%>%
  ungroup()

cortex_d <- diameter%>%
  filter(type == "cortex")
endo_d <- diameter%>%
  filter(type == "endodermis")
stele_d <- diameter%>%
  filter(type == "stele")

diam <- rbind(stele_d, cortex_d)%>%
  mutate(d_cortex = c(rep(NA, nrow(stele_d)),c(cortex_d$d-endo_d$d)),
         tissue = type)

Data_params%>%
  filter(tissue %in% c("stele", "cortex"))%>%
  ggplot(aes(dist, layer))+
  geom_point(aes(colour = Root), alpha = 0.1)+
  geom_smooth(aes(colour = Root), method = "lm")+
  geom_point(aes(dist, d*1000, colour = name), data = diam%>%
               filter(tissue == "stele"), shape = 2)+
  geom_point(aes(dist, d_cortex*500, colour = name), data = diam%>%
               filter(tissue == "cortex"), shape = 2)+
  facet_wrap(~tissue)+
  labs(shape = "GRANAR")+
  viridis::scale_colour_viridis(discrete= T)

MECHA

From the generated anatomies the MECHA model estimates the kr with cell hydraulic properties set in the input file.

fls <- list.files("./MECHA_GRANAR/cellsetdata/")
fls <- fls[grepl(".xml", fls)]
fls <- fls[fls != "current_root.xml"]

# fls <- fls[fls %!in% fls_done]

fls_done <- NULL

for(j in fls){
  print(j)
  id_root <- unlist(str_split(j, "_"))[1]
  id_dist <- unlist(str_split(j, "_"))[2]
  id_rep <- parse_number(unlist(str_split(j, "_"))[3])
  if(id_dist == "A"){
    id_root <- "LAT_A"
    id_dist <- unlist(str_split(j, "_"))[3]
    id_rep <- parse_number(unlist(str_split(j, "_"))[4])
  }
  
  if(file.exists("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Project_Test/results/Macro_prop_1,0.txt")){
    file.remove("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Project_Test/results/Macro_prop_1,0.txt")
    file.remove("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Project_Test/results/Macro_prop_2,1.txt")
    file.remove("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Project_Test/results/Macro_prop_4,2.txt")
    file.remove("./MECHA_GRANAR/cellsetdata/current_root.xml")
  }
  
  fc <- file.copy(paste0("./MECHA_GRANAR/cellsetdata/",j), "./MECHA_GRANAR/cellsetdata/current_root.xml", overwrite = T)
  if(fc){
  system("C:/Users/heymansad/AppData/Local/Continuum/anaconda3/envs/MECHA/python.exe MECHA_GRANAR/MECHAv4_GRANAR.py")
  if(file.exists("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Project_Test/results/Macro_prop_1,0.txt")){
    file.copy("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Project_Test/results/Macro_prop_1,0.txt",
              paste0("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Macro_prop_1,0_",id_root,"_",id_dist,"_",id_rep,".txt"))
    file.copy("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Project_Test/results/Macro_prop_2,1.txt",
              paste0("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Macro_prop_2,1_",id_root,"_",id_dist,"_",id_rep,".txt"))
    file.copy("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Project_Test/results/Macro_prop_4,2.txt",
              paste0("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Macro_prop_4,2_",id_root,"_",id_dist,"_",id_rep,".txt"))
    tmp <- paste0(id_root, "_", id_dist, "_", id_rep,".xml")
    fls_done <- c(fls_done, tmp)
  }else{
    print(j)
    print("fail")
  }
  }else{break}
  
}

Loading the output of MECHA.

K <- NULL

for(h in fls_done){
  
  id_root <- unlist(str_split(h, "_"))[1]
  id_dist <- unlist(str_split(h, "_"))[2]
  id_rep <- parse_number(unlist(str_split(h, "_"))[3])
  if(id_dist == "A"){
    id_root <- "LAT_A"
    id_dist <- unlist(str_split(h, "_"))[3]
    id_rep <- parse_number(unlist(str_split(h, "_"))[4])
  }

  path1 <- paste0("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Macro_prop_1,0_",id_root,"_",id_dist,"_",id_rep,".txt")
  path2 <- paste0("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Macro_prop_2,1_",id_root,"_",id_dist,"_",id_rep,".txt")
  path3 <- paste0("./MECHA_GRANAR/Projects/GRANAR/out/M1v4/Root/Macro_prop_4,2_",id_root,"_",id_dist,"_",id_rep,".txt")
  
  tmp_K <- read_mecha_output(path1, path2, path3)
  tmp_K <- cbind(root_name = id_root, d = parse_number(id_dist), rep = id_rep, tmp_K)
  K <- rbind(K, tmp_K)
}

Comparison between the hydraulic conductivity curve that were generated and the ones of Doussan et al. 1998.

K_xyl_unM <- nodes%>%
  filter(type == "xylem")%>%
  dplyr::group_by(id_cell, repet, name, dist)%>%
  dplyr::summarise(a = mean(area)*1000^2, ## mm2 to um2
            k_xyl_s = a^2/(8*pi*200*1E-5/3600/24)*1E-12)%>%
  ungroup()%>%
  dplyr::group_by(repet, name, dist)%>%
  dplyr::summarise(Kx_un = sum(k_xyl_s)*200/1E4)%>%
  ungroup()%>%
  mutate(ID = paste0(name,"_",dist,"_", repet))
K$ID = paste0(K$root_name,"_",K$d, "_", K$rep)

Kconduct <- left_join(K, K_xyl_unM, by = "ID")

# Load Doussan et al. 1998 Conductivities data file
conductivities <- read.csv("Doussan_conductivities.csv")
conv = 0.001157407
conv_kx = 1.157407e-09

K_axial <- NULL
for(root in c("SBR", "Tap", "Basal", "LAT", "LAT_A")){
  if (root == "SBR"){
    unM = 10
    M = 12
  }else if (root == "Basal"){
    unM = 15
    M = 20
  }else if(root == "Tap"){
    unM = 10
    M = 12
  }else if (root == "LAT"){
    unM = 1
    M = 2
  }else if (root == "LAT_A"){
    unM = 5
    M = 7
  }
  
  kx1 <- tibble(root_name = root, d = Kconduct$d[Kconduct$d <= unM & Kconduct$root_name == root], k = Kconduct$Kx_un[Kconduct$d <= unM & Kconduct$root_name == root])
  kx2 <- tibble(root_name = root, d = Kconduct$d[Kconduct$d >= M & Kconduct$root_name == root], 
                  k = Kconduct$kx[Kconduct$d >= M & Kconduct$root_name == root])
  
  KX <- rbind(kx1,kx2)
  K_axial <- rbind(K_axial, KX)
}



K_axial %>%
  ggplot()+
  geom_point(aes(d, k*conv_kx, colour = root_name))+
  geom_line(aes(d, k*conv_kx, colour = root_name), size = 1, alpha = 0.7)+
  #geom_line(aes(dist,K, colour = Root), size = 1.5, alpha = 0.5, data = conductivities%>%filter(type == "kx"))+
  ylab("Axial conductance [m4 s-1 MPa-1]")+
  xlim(0,45)+
  scale_y_continuous(limits =  c(0.8E-14,1E-8), trans = 'log10')+
  viridis::scale_colour_viridis(discrete= T)

K_axial %>%
  ggplot()+
  #geom_point(aes(d, k*conv_kx, colour = root_name))+
  geom_line(aes(d, k*conv_kx, group = root_name), size = 1, alpha = 0.3)+
  geom_line(aes(dist,K, linetype = Root), size = 1.5, alpha = 1, data = conductivities%>%filter(type == "kx"))+
  ylab("Axial conductance [m4 s-1 MPa-1]")+
  xlim(0,100)+
  scale_y_continuous(trans = 'log10')+
  viridis::scale_colour_viridis(discrete= T)

K_radial <- NULL
for(root in c("SBR", "Tap", "Basal", "LAT", "LAT_A")){
  if (root == "SBR"){
    endo_sub = 12
    exo_casp = 18
  }else if (root == "Basal"){
    endo_sub = 10
    exo_casp = 15
  }else if(root == "Tap"){
    endo_sub = 8
    exo_casp = 13
  }else if (root == "LAT"){
    endo_sub = 3
    exo_casp = 5
  }else if (root == "LAT_A"){
    endo_sub = 4
    exo_casp = 6
  }
  
  kr1 <- tibble(root_name = root, d = K$d[K$d <= endo_sub & K$root_name == root], k = K$kr1[K$d <= endo_sub & K$root_name == root])
  kr2 <- tibble(root_name = root, d = K$d[K$d <= exo_casp & K$d > endo_sub & 
                                                  K$root_name == root ], 
                  k = K$kr2[K$d <= exo_casp & K$d > endo_sub & K$root_name == root])
  kr3 <- tibble(root_name = root, d = K$d[K$d > exo_casp & K$root_name == root ], k = K$kr3[K$d > exo_casp & K$root_name == root])
  
  KR <- rbind(kr1,kr2, kr3)
  K_radial <- rbind(K_radial, KR)
}



K_radial %>%
  ggplot()+
  geom_point(aes(d, k*conv, colour = root_name))+
  geom_line(aes(d, k*conv, colour = root_name), size = 1)+
  #geom_line(aes(dist,K, colour = Root), size = 1.5, alpha = 1, data = conductivities%>%filter(type == "kr"))+
  ylab("Radial conductivity [m s-1 MPa-1]")+
  xlim(0,45)+
  viridis::scale_colour_viridis(discrete= T)

K_radial %>%
  ggplot()+
  #geom_point(aes(d, k*conv), alpha = 0.3)+
  geom_line(aes(d, k*conv, group = root_name), size = 1, colour = "black", alpha = 0.3)+
  geom_line(aes(dist,K, linetype = Root), size = 1.5, alpha = 1, data = conductivities%>%filter(type == "kr"))+
  ylab("Radial conductivity [m s-1 MPa-1]")+
  xlim(0,100)+
  viridis::scale_colour_viridis(discrete= T)